---
title: |
  Replication material: Anghel and Schulte-Cloos (2022) 'COVID-19 related anxieties do not decrease support for liberal democracy, Harvard Dataverse DOI: https://doi.org/10.7910/DVN/QKLSWY.
author: |
  Julia Schulte-Cloos
date: |
  `r lubridate::today()`
output:
  bookdown::html_document2:
    theme: flatly
    highlight: textmate
    toc: yes
    toc_float:
      collapsed: no
      smooth_scroll: no
    number_sections: yes
    code_folding: hide
    toc_depth: 2 
link-citations: true
---

```{r setup}

pacman::p_load(
  tidyverse,
  tidymodels,
  rsample,
  knitr,
  kableExtra,
  RColorBrewer,
  tictoc,
  gridExtra,
  Hmisc, 
  colorspace,
  scales,
  broom,
  jcolors,
  lubridate,
  showtext,
  modelsummary,
  tidyr,
  embed,
  ggdist,
  distributional, 
  xfun, 
  benchmarkme
)


# set seed specified in the PAP
set.seed(853147)

knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  cache = FALSE,
  echo = TRUE, 
  fig.retina = 2
)


# evaluate computationally intense code chunks? 
eval_computationally_intense <- FALSE

```



```{r set-ggplot-font-html}

# add a font from google fonts
font_add_google(
  name = "Fira Sans",
  family = "Fira Sans"
)

ggplot_font <- "Fira Sans"
```


```{r theme-set-ggplot}

ggplot2::theme_set(
  theme_bw(base_size = 24) +
    theme(text = element_text(family = ggplot_font)) +
    theme(
      plot.title = element_text(
        hjust = 0.5
        ),
      plot.subtitle = element_text(
        size = rel(0.8), 
        hjust = 0.5
        ), 
      legend.title = element_text(
        size = rel(0.8)
        ),
      plot.caption = element_text(hjust = 0, 
                              size = rel(0.7))
      ))


showtext::showtext_auto()
```


# Data Import and Operationalization 

## Survey experimental data

We import the (raw) survey experimental data. 

```{r load-raw-data}

covid_fear_data <- read_csv("data/raw_covid_fear.csv")
```


## Recoding of covariates

We recode some of the variables. 

```{r recode-covariates}

covid_fear_data <-
  covid_fear_data %>%
  mutate(
    gender = case_when(
      gender == 1 ~ 1,
      gender == 2 ~ 0,
      gender == 3 ~ NA_real_
    ),
    age = as.numeric(2021 - year_born),
    diaspora = case_when(
      diaspora == 1 ~ 0,
      diaspora == 2 ~ 1,
      diaspora == 3 ~ NA_real_
    ),
    minority = case_when(
      minority == 1 ~ 1,
      minority == 2 ~ 0,
      minority == 3 ~ NA_real_
    ),
    religion_life = case_when(
      religion_1 == 1 ~ 1,
      religion_1 == 2 ~ 0,
      religion_1 == 3 ~ NA_real_
    ),
    religious_service = case_when(
      religion_2 == 6 ~ NA_real_,
      TRUE ~ as.numeric(religion_2)
    ),
    covid_infection = case_when(
      covid_infection == 1 ~ 1,
      covid_infection == 2 ~ 0
    ),
    covid_vaccination = case_when(
      covid_vaccination == 1 ~ 1,
      covid_vaccination == 2 ~ 0
    )
  ) %>%
  rename(female = gender) %>%
  select(-religion_1, -religion_2)
```


## Aggregate statistical data

We also add data on the official numbers of COVID-19 infections during our field work (weekly data obtained by the ECDC).

```{r add-ecdc-covid-infections-data}

covid_fear_data <- covid_fear_data %>%
  left_join(.,
    read_csv("./data/covid_incidence_ecdc_data.csv") %>%
      mutate(nuts2_code = str_sub(nuts_code, start = 1, end = 4)) %>%
      group_by(nuts2_code, year_week) %>%
      summarise(rate_14_day_per_100k = mean(
        rate_14_day_per_100k,
        na.rm = T
      )) %>%
      rename(response_week_incidence = rate_14_day_per_100k),
    by = c(
      "region_nuts" = "nuts2_code",
      "survey_year_week" = "year_week"
    )
  )
```


## Dimensionality reduction

We first subset the data to include only the relevant outcome variables, the covariates, the treatment indicator, and the relevant identifiers. We then perform the dimensionality reduction of the different outcome variables (by country). We reduce the items into a single index by performing a Principal Component Analysis (PCA), by relying on an supervised method (UMAP), and by averaging across all items (Means). 

```{r outcome-covariates-subset-data}

outcome_variables <- covid_fear_data %>%
  select(starts_with("rwa_") |
    starts_with("outgroup_") |
    starts_with("nationalism_") |
    starts_with("covid_policy_authoritarianism_") |
    starts_with("covid_policy_nationalism_") |
    starts_with("covid_policy_outgroup_")) %>%
  colnames()


covariates <- c(
  "female",
  "age",
  "diaspora",
  "minority",
  "religion_life",
  "covid_infection",
  "covid_vaccination",
  "education_level",
  "urban_rural",
  "religious_service",
  "support_government",
  "political_news",
  "response_week_incidence"
)


covid_fear_data <- covid_fear_data %>%
  ungroup() %>%
  select(
    "ID",
    "country",
    "region_nuts",
    "fear_treatment",
    all_of(covariates),
    all_of(outcome_variables),
    starts_with("emotion")
  )
```

Next, we nest the data by the different outcome dimensions.

```{r nested-by-outcome-dimensions}

covid_dimensions <- covid_fear_data %>%
  ungroup() %>%
  select(
    country, ID,
    matches("covid_policy_authoritarianism_"),
    matches("covid_policy_outgroup_"),
    matches("covid_policy_nationalism_"),
    matches("rwa_"),
    matches("outgroup_"),
    matches("nationalism_")
  ) %>%
  pivot_longer(starts_with(c(
    "rwa_", "outgroup_",
    "nationalism_", "covid_policy_"
  )),
  names_to = "outcome_variable",
  values_to = "value"
  ) %>%
  mutate(outcome_group = case_when(
    str_detect(outcome_variable, "covid_policy_outgroup_") ~ "outgroup_covid",
    str_detect(outcome_variable, "covid_policy_nationalism_") ~ "nationalism_covid",
    str_detect(outcome_variable, "covid_policy_authoritarianism") ~ "rwa_covid",
    str_detect(outcome_variable, "rwa_") ~ "rwa",
    str_detect(outcome_variable, "outgroup_") ~ "outgroup",
    str_detect(outcome_variable, "nationalism_") ~ "nationalism"
  )) %>%
  group_by(country, outcome_group) %>%
  nest() %>%
  # bring back to wide format for the dimensionality reduction
  mutate(data = map(
    data,
    ~ pivot_wider(.,
      id_cols = "ID",
      names_from = "outcome_variable",
      values_from = "value"
    )
  )) %>%
  # calculate mean of each outcome group
  mutate(mean_group = map(
    data,
    ~ .x %>%
      select(-ID) %>%
      rowMeans()
  ))
```


### PCA

We perform PCA on each of the dimensions. 

```{r pca-dimensionality-reduction}

covid_pca <- covid_dimensions %>%
  # define recipe, identifier ID, normalise, and apply umap algorithm
  mutate(
    pca_rec = map(
      data,
      ~ recipe(~.,
        data = .x
      ) %>%
        update_role(ID, new_role = "id") %>%
        step_normalize(all_predictors()) %>%
        step_pca(all_predictors(),
          options = c(
            center = TRUE,
            scale. = TRUE
          )
        )
    ),
    # preprocess pca
    pca_prep = map(
      pca_rec,
      ~ prep(.)
    ),
    # add pca to the data
    pca_bake = map2(
      pca_prep,
      data,
      ~ bake(.x,
        new_data = .y
      )
    )
  )


# direction of pca
covid_pca <- covid_pca %>%
  mutate(pca_1 = map(
    pca_bake,
    ~ pull(., PC1)
  )) %>%
  mutate(transformation = map2(
    pca_1,
    mean_group,
    ~ ifelse(
      cor(.x, .y, method = c("pearson")) < 0,
      TRUE,
      FALSE
    )
  )) %>%
  mutate(pca_first = case_when(
    transformation == TRUE ~ map(pca_1, ~ .x * (-1)),
    transformation == FALSE ~ pca_1
  )) %>%
  mutate(correlation = map2(
    pca_first,
    mean_group,
    ~ cor(.x, .y, method = c("pearson"))
  ))


# create pca_wide dataset
pca_first_wide <- covid_pca %>%
  # add the ID value
  mutate(data = map(
    data,
    ~ select(., ID)
  )) %>%
  select(outcome_group, pca_first, data) %>%
  unnest(pca_first, data) %>%
  pivot_wider(
    id_cols = c("ID"),
    names_from = "outcome_group",
    names_prefix = "pca_",
    values_from = "pca_first"
  )


# add pca to the data
covid_fear_data <- covid_fear_data %>%
  left_join(
    .,
    pca_first_wide
  ) %>%
  group_by(country) %>%
  mutate(across(starts_with("pca_"), ~ as.numeric(scale(.x))))


```

Table \@ref(tab:pca-statistics) shows the metrics of the PCA that we conducted to arrive at the conceptually relevant outcome dimensions of interest. We show the first two components and the amounts of variance in the initial set of variables that can be explained by these three components. The percentage of variance explained by each principal component is the eigenvalue of each component divided by the sum of all eigenvalues.

```{r pca-statistics}

pca_stats <- covid_pca %>%
  select(pca_prep) %>%
  mutate(prcomp = map(
    pca_prep,
    # get pr comp object
    ~ pluck(., "steps") %>%
      .[[2]] %>%
      pluck("res")
  )) %>%
  mutate(tidy_prcomp = map(
    prcomp,
    ~ tidy(.,
      matrix = "eigenvalues"
    ) %>%
      mutate(N_items = n()) %>%
      slice(1:2)
  )) %>%
  arrange(outcome_group, country) %>%
  unnest(tidy_prcomp) %>%
  arrange(desc(outcome_group)) %>%
  mutate(
    type = case_when(
      str_detect(
        outcome_group,
        "covid"
      ) ~ "COVID-19 Policies",
      TRUE ~ "General Attitudes"
    ),
    dimension = str_to_title(str_remove_all(
      outcome_group,
      "pca_|_covid"
    )),
    dimension = case_when(
      dimension == "Rwa" ~ "Authoritarian",
      dimension == "Nationalism" ~ "Nationalist",
      dimension == "Outgroup" ~ "Outgroup-Hostile",
      TRUE ~ dimension
    )
  ) %>%
  arrange(type, desc(outcome_group)) %>%
  select(-pca_prep, -prcomp)


pca_stats %>%
  ungroup() %>%
  select(country, dimension, PC, N_items, percent) %>%
  # variance explained by each PC
  kable(.,
    col.names = c("Country", "Dimension", "Principal Component", "N Items", "% Variance Explained"),
    digits = 2,
    booktabs = T,
    caption = "First and second principal components of each conceptually relevant dimension and amount of variance explained by each component."
  ) %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover"),
    font_size = 11
  ) %>%
  pack_rows(index = table(pca_stats$type)) %>%
  row_spec(which(pca_stats$PC == 1), bold = T)
```


### UMAP

We rely on the UMAP algorithm for another way of reducing the dimensionality in the data.

```{r umap-dimensionality-reduction, eval=eval_computationally_intense}

# for computational efficiency, we do not call the UMAP algorithm each time, but store the data in a csv, which we add to the other data in the following code chunk.

set.seed(853147)

covid_umap <- covid_dimensions %>%
  # define recipe, identifier ID, normalise, and apply umap algorithm
  mutate(
    umap_rec = map(
      data,
      ~ recipe(~.,
        data = .x
      ) %>%
        update_role(ID, new_role = "id") %>%
        step_normalize(all_predictors()) %>%
        step_umap(all_predictors(),
          seed = c(
            853147,
            853147
          ),
          num_comp = 1
        )
    ),
    # preprocess umap
    umap_prep = map(
      umap_rec,
      ~ prep(.)
    ),
    # add umap to the data
    umap_bake = map2(
      umap_prep,
      data,
      ~ bake(.x,
        new_data = .y
      )
    )
  )


# direction of umap
covid_umap <- covid_umap %>%
  mutate(umap_1 = map(
    umap_bake,
    ~ pull(., -ID)
  )) %>%
  mutate(transformation = map2(
    umap_1,
    mean_group,
    ~ ifelse(
      cor(.x, .y, method = c("pearson")) < 0,
      TRUE,
      FALSE
    )
  )) %>%
  mutate(umap_first = case_when(
    transformation == TRUE ~ map(umap_1, ~ .x * (-1)),
    transformation == FALSE ~ umap_1
  )) %>%
  mutate(correlation = map2(
    umap_first,
    mean_group,
    ~ cor(.x, .y, method = c("pearson"))
  ))


# create umap_wide dataset
umap_first_wide <- covid_umap %>%
  # add the ID value
  mutate(data = map(
    data,
    ~ select(., ID)
  )) %>%
  select(outcome_group, umap_first, data) %>%
  unnest(umap_first, data) %>%
  pivot_wider(
    id_cols = c("ID"),
    names_from = "outcome_group",
    names_prefix = "umap_",
    values_from = "umap_first"
  )

# save umap data for computational efficiency
umap_first_wide %>%
  write_csv("./data/covid_umap.csv")
```

```{r add-umap-standardise}

# add umap to the data
covid_fear_data <- covid_fear_data %>%
  left_join(
    .,
    read_csv("./data/covid_umap.csv")
  ) %>%
  group_by(country) %>%
  mutate(across(starts_with("umap_"), ~ as.numeric(scale(.x))))
```


### Means 

Finally, we also compute the simple means of each items that belong to the different conceptually relevant outcome dimensions of interest.

```{r means-dimensionality-reduction}

covid_means <- covid_dimensions %>%
  select(outcome_group, data) %>%
  mutate(mean = map(
    data,
    ~ select(., -ID) %>%
      rowMeans(.)
  ))

means_wide <- covid_means %>%
  mutate(data = map(
    data,
    ~ select(., ID)
  )) %>%
  unnest(mean, data) %>%
  pivot_wider(
    id_cols = c("country", "ID"),
    names_from = "outcome_group",
    names_prefix = "mean_",
    values_from = "mean"
  )


# add means to the data
covid_fear_data <- covid_fear_data %>%
  left_join(
    .,
    means_wide
  ) %>%
  group_by(country) %>%
  mutate(across(starts_with("mean_"), ~ as.numeric(scale(.x))))
```


# Data Analysis 

## Treatment effectiveness

We first assess the effectiveness of the treatment by analysing the results of our manipulation check included *after* respondents answer all the questions related to our outcome variables of interest. 

### Figure 1: Emotional states between Treatment and Control

Figure \@ref(fig:manipulation-check) shows respondents' emotions when thinking about infectious diseases like COVID-19 by their treatment status. Table \@ref(tab:means-standard-deviations-emotions-unscaled) shows the group means and respective differences in the unstandardised data.

```{r manipulation-check, fig.cap="Means of emotional responses among treatment and control group when thinking about infectious diseases like COVID-19."}

# levels of emotionality among treated and non-treated: long data (standardised)
manipulation_check <- covid_fear_data %>%
  ungroup() %>%
  select(
    fear_treatment, country,
    starts_with("emotions")
  ) %>%
  # standardise with mean 0 and sd 1
  mutate(across(
    .cols = -c("fear_treatment", "country"),
    scale
  )) %>%
  pivot_longer(starts_with("emotions"),
    names_to = "emotion",
    values_to = "value",
    names_prefix = "emotions_"
  ) %>%
  mutate(emotion = str_to_title(emotion))



manipulation_check_plot <- manipulation_check %>%
  left_join(.,
    # left_join the augmented data from fitting each regression
    manipulation_check %>%
      group_by(emotion) %>%
      nest() %>%
      mutate(
        lm_emotion =
          map(
            data,
            ~ lm(value ~ fear_treatment, data = .)
          ),
        fitted_emotion =
          map(
            lm_emotion,
            ~ broom::augment(.x,
              se_fit = TRUE,
              interval = c("confidence")
            )
          )
      ) %>%
      unnest(fitted_emotion),
    by = c("fear_treatment", "emotion")
  ) %>%
  distinct(fear_treatment, emotion, .fitted, .se.fit, .lower, .upper)


# estimate one model for the extraction of degrees of freedom
m_emotions <- lm(value ~ fear_treatment, data = manipulation_check %>%
  filter(emotion == "Hope"))

# calculate the positions of the arrows and labels that we use to highlight two emotions
overallmean_treatment <- manipulation_check_plot %>%
  filter(emotion == "Worried" & fear_treatment == 1) %>%
  pull(.fitted)
overallmean_control <- manipulation_check_plot %>%
  filter(emotion == "Afraid" & fear_treatment == 0) %>%
  pull(.fitted)


# coordinates for arrows
arrows <- tibble(
  xfinish = c(0.15, -0.15),
  fear_treatment = c(1, 0),
  xarrow = c(
    overallmean_treatment + 0.03,
    overallmean_control - 0.03
  ),
  yarrow = c(6, 1.2),
  yfinish = c(5.2, 2.2)
)

arrows <- arrows %>%
  mutate(
    ytext = ifelse(fear_treatment == 1, yfinish - 0.2,
      yfinish + 0.2
    ),
    xtext = xfinish
  )

scale_shape_values <- c("Treatment" = 10,
                        "Control" = 8)

# render the final plot
manipulation_check_plot %>%
  mutate(fear_treatment = ifelse(fear_treatment == 1,
    "Treatment",
    "Control"
  )) %>%
  ggplot(aes(y = emotion)) +
  stat_dist_halfeye(
    aes(
      dist = dist_student_t(
        df = df.residual(m_emotions),
        mu = .fitted,
        sigma = .se.fit
      ),
      color = as.factor(fear_treatment), 
      shape = as.factor(fear_treatment)
    ),
    size = 6,
    fill = "gray95",
    position = position_dodge(width = 0.5)
  ) +
  # add arrows
  geom_curve(
    data = arrows,
    aes(
      x = xfinish,
      xend = xarrow,
      y = yfinish,
      yend = yarrow
    ),
    arrow = arrow(length = unit(0.08, "inch")),
    color = "gray80",
    curvature = 0.3
  ) +
  annotate("text",
    y = arrows$ytext,
    x = arrows$xtext,
    label = arrows %>%
      mutate(fear_treatment = ifelse(fear_treatment == 1,
        "Treatment",
        "Control"
      )) %>%
      pull(fear_treatment),
    size = 8,
    lineheight = .2, 
    family = ggplot_font
  ) +
  xlab("") +
  ylab("") +
  scale_x_continuous(limits = c(-0.2, 0.2)) +
  labs(subtitle = "When you think about infectious diseases like COVID-19, to what extent do you\nexperience the following emotions?") +
  jcolors::scale_color_jcolors(
    name = "",
    palette = "pal9",
    guide = "none"
  ) + 
  scale_shape_manual(values = scale_shape_values, 
                     guide = "none")
```


### Heat Map of Emotion States

Figure \@ref(fig:corr-emotions) shows a heat map of the bivariate correlations between the different emotional states that respondents reported in the treatment and control group. We find strong positive correlations between the emotional states of feeling worried and afraid and between feeling angry and outraged (pearson's r > 0.5). There is a modest positive correlation between feeling angry and worried (pearson's r > 0.3) and a modest negative correlation between between feeling hopeful and afraid (pearson's r < -0.3). There are no statistically significant differences in these general patterns among treatment and control group. 

```{r corr-emotions, fig.width=10, fig.cap="Heat map of correlations between different emotional states among treated and control respondents."}

# correlation matrix by treatment/control
emotion_correlation <- bind_rows(
  cor(covid_fear_data %>%
    ungroup() %>%
    filter(fear_treatment == 1) %>%
    select(starts_with("emotion"))) %>%
    as_tibble(.) %>%
    mutate(emotion_1 = colnames(.)) %>%
    # bring into long format
    pivot_longer(
      cols = -c("emotion_1"),
      names_to = "emotion_2",
      values_to = "value"
    ) %>%
    mutate(emotion_1 = str_to_title(
      str_remove(
        emotion_1,
        "emotions_"
      )
    )) %>%
    mutate(emotion_2 = str_to_title(
      str_remove(
        emotion_2,
        "emotions_"
      )
    )) %>%
    mutate(condition = "Fear Condition"),
  cor(covid_fear_data %>%
    ungroup() %>%
    filter(fear_treatment == 0) %>%
    select(starts_with("emotion"))) %>%
    as_tibble(.) %>%
    mutate(emotion_1 = colnames(.)) %>%
    # bring into long format
    pivot_longer(
      cols = -c("emotion_1"),
      names_to = "emotion_2",
      values_to = "value"
    ) %>%
    mutate(emotion_1 = str_to_title(
      str_remove(
        emotion_1,
        "emotions_"
      )
    )) %>%
    mutate(emotion_2 = str_to_title(
      str_remove(
        emotion_2,
        "emotions_"
      )
    )) %>%
    mutate(condition = "Control Condition")
)



ggplot(
  data = emotion_correlation %>%
    filter(!is.na(value)) %>%
    mutate(value = round(value, digits = 2)),
  aes(
    x = emotion_2,
    y = fct_rev(emotion_1),
    fill = value
  )
) +
  geom_tile(
    color = "white",
    lwd = 1.5,
    linetype = 1
  ) +
  geom_text(aes(label = value),
    color = "white",
    size = 9
  ) +
  labs(
    x = "",
    y = ""
  ) +
  scale_fill_jcolors_contin() +
  theme(legend.position = "none") +
  facet_wrap(~ fct_rev(condition))
```


### Group Differences in Emotion States 
 
In addition to the visual presentation of group differences in respondents' (standardized) manipulated emotions presented in Figure \@ref(fig:manipulation-check), Table \@ref(tab:means-standard-deviations-emotions-unscaled) shows the group means and respective differences in the unstandardized data.

```{r means-standard-deviations-emotions-unscaled}


manipulation_check_unscaled <- covid_fear_data %>%
  select(
    ID, fear_treatment,
    starts_with("emotions")
  ) %>%
  pivot_longer(starts_with("emotions"),
    names_to = "emotion",
    values_to = "value",
    names_prefix = "emotions_"
  ) %>%
  mutate(emotion = str_to_title(emotion)) %>%
  pivot_wider(
    names_from = emotion,
    values_from = value,
    id_cols = c(ID, fear_treatment)
  ) %>%
  mutate(fear_treatment = if_else(fear_treatment == 0, "Control", "Treatment"))

caption <- "Means and standard deviations of manipulated emotions among treated and control respondents"


datasummary_balance(~fear_treatment,
  data = manipulation_check_unscaled,
  caption = caption,
  fmt = 2,
  output = "kableExtra"
)
```



## Treatment effects estimation

Next, we estimate the effects of fear of COVID-19 on higher-level attitudes and COVID-19 specific preferences. 
To do so, we bring the data in a long format which we then nest by outcome dimensions (`covid_fear_nested` and also by countries `covid_fear_nested_by_countries`) for a more efficient estimation of the different models.

```{r data-nested-by-outcome-dimensions}

categorical_covariates <- c(
  "female",
  "diaspora",
  "minority",
  "religion_life",
  "covid_infection",
  "covid_vaccination",
  "region_nuts"
)
numerical_covariates <- c(
  "age",
  "education_level",
  "urban_rural",
  "religious_service",
  "support_government",
  "political_news",
  "response_week_incidence"
)


# pivot longer of pca, means, and umap
covid_fear_nested <- covid_fear_data %>%
  # change the names of the covid dimensionality reduction for an easy pivot_longer
  rename_with(
    cols = matches("_covid"),
    .fn = ~ str_replace_all(.,
      pattern = "_covid",
      replacement = "-covid"
    )
  ) %>%
  ungroup() %>%
  select(
    "ID",
    "fear_treatment",
    "country",
    all_of(categorical_covariates),
    all_of(numerical_covariates),
    starts_with("pca_"), starts_with("umap_"), starts_with("mean_")
  ) %>%
  pivot_longer(
    cols = c(
      starts_with("pca_"),
      starts_with("umap_"),
      starts_with("mean")
    ),
    names_to = c(".value", "variable"),
    names_pattern = "(.*)_(.*)"
  ) %>%
  mutate(variable = str_replace_all(variable,
    pattern = "-",
    replacement = "_"
  )) %>%
  rename_with(
    .cols = c("pca", "mean", "umap"),
    .fn = ~ paste0(., "_value")
  )



# keep a separate data that is nested by country
covid_fear_nested_by_country <- covid_fear_nested %>%
  ungroup() %>%
  nest(data = c(-variable, -country)) %>%
  mutate(data = map(
    data,
    ~ mutate(
      .,
      across(
        all_of(categorical_covariates),
        ~ as.factor(.)
      )
    )
  ))

# nest the data by outcome dimensions
covid_fear_nested <- covid_fear_nested %>%
  nest(data = c(-variable)) %>%
  mutate(data = map(
    data,
    ~ mutate(
      .,
      across(
        all_of(categorical_covariates),
        ~ as.factor(.)
      )
    )
  ))
```

We define the linear regression model to be estimated on each model as a function that takes five arguments. These arguments specify whether we want to estimate a simple or a full model (`model`), what the value to be predicted should be (i.e., the outcome dimension, `outcome`), and the data to be used for estimation (`data`). Finally, the function takes a logical argument indicating whether we will directly tidy the model output (`tidy`) and another one indicating whether the model should be estimated without country fixed effects, which we do not need when estimating the model separately by each country (`by_country`).

```{r define-regression-function}

# lm function with country-fixed effects
fit_lm <- function(model,
                   outcome,
                   data,
                   tidy = FALSE,
                   by_country = FALSE) {
  data <- data
  model <- model
  outcome <- outcome
  covariates_in_data <- covariates

  # subset data (depending on model - simple vs. fully specified)
  if (model == "simple") {
    data <- data %>%
      select(any_of(c(
        regression_outcome_variables, "fear_treatment",
        "country", "region_nuts"
      )))
  } else if (model == "full") {
    data <- data %>%
      select(any_of(c(
        regression_outcome_variables, "fear_treatment", 
        "country", "region_nuts",
        covariates
      ))) %>%
      # omit factors in the bootstrap samples that contain just one level
      janitor::remove_constant(na.rm=T)
    # vector with all covariates
    covariates_in_data <- colnames(data %>% select(any_of(covariates)))
  }


  # lm simple model
  if (by_country == FALSE) {
    lm_formula_simple <- as.formula(
      paste(outcome,
        paste0(c("fear_treatment", "country"), collapse = " + "),
        sep = " ~ "
      )
    )
  } else if (by_country == TRUE) {
    lm_formula_simple <- as.formula(
      paste(outcome,
        paste0(c("fear_treatment"), collapse = " + "),
        sep = " ~ "
      )
    )
  }

  # lm full model
  if (by_country == FALSE) {
    lm_formula_full <- as.formula(
      paste(outcome,
        paste0(c(
          "fear_treatment",
          covariates_in_data,
          "country"
        ), collapse = " + "),
        sep = " ~ "
      )
    )
  } else if (by_country == TRUE) {
    lm_formula_full <- as.formula(
      paste(outcome,
        paste0(c(
          "fear_treatment",
          covariates_in_data
        ), collapse = " + "),
        sep = " ~ "
      )
    )
  }


  # lm formula (depending on model - simple vs. fully specified)
  lm_formula <- if (model == "simple") {
    lm_formula_simple
  } else if (model == "full") {
    lm_formula_full
  }

  # estimate model
  if (tidy == FALSE) {
    return(
      eval(bquote(lm(.(lm_formula), data = data)))
    )
  } else if (tidy == TRUE) {
    return(
      eval(bquote(lm(.(lm_formula), data = data))) %>%
        tidy(conf.int = TRUE)
    )
  }
}
```

```{r regression-variables}

regression_outcome_variables <- c(
  "pca_value",
  "umap_value",
  "mean_value"
)

covariates <- c(
  "female",
  "age",
  "diaspora",
  "minority",
  "religion_life",
  "covid_infection",
  "covid_vaccination",
  "education_level",
  "urban_rural",
  "religious_service",
  "support_government",
  "political_news",
  "response_week_incidence"
)
```

#### Bootstrap estimation 

Next, we estimate the bootstrapped confidence intervals by repeatedly resampling from our data (5000 times) and re-estimating the model on each of the splits. 

```{r estimation-bootstrap-confidence-intervals, eval=eval_computationally_intense}

set.seed(853147)
# for computational efficiency, we save the bootstrapped CIs in a csv
boot_times <- 5000

tic()
bootstraps <- covid_fear_nested %>%
  # sample 5000 times over each subset of the data
  mutate(boots = map(
    .x = data,
    .f = ~ rsample::bootstraps(.x,
      strata = "country",
      times = boot_times,
      apparent = T
    )
  )) %>%
  # estimate model on each of the splits
  mutate(boots = map(
    .x = boots,
    ~ mutate(
      .x,
      # estimate lm model simple
      lm_model_simple_tidy = pmap(
        list(..1 = .x$splits),
        .f = ~ map(
          regression_outcome_variables,
          function(regression_outcome_variables) {
            fit_lm(
              model = "simple",
              outcome = paste(regression_outcome_variables),
              data = analysis(..1),
              tidy = TRUE
            )
          }
        )
      ),
      # add names
      lm_model_simple_tidy = map(
        lm_model_simple_tidy,
        ~ set_names(.,
          nm = c("PCA", "UMAP", "MEAN")
        )
      ),
      # turn list of tibbles into one tibble and create unique names in the term column required for int_pctl
      lm_model_simple_tidy = map(
        lm_model_simple_tidy,
        ~ map_dfr(.,
          bind_rows,
          .id = "regression_outcome"
        ) %>%
          mutate(term = paste(regression_outcome, term, sep = "."))
      ),
      # estimate lm model full
      lm_model_full_tidy = pmap(
        list(..1 = .x$splits),
        .f = ~ map(
          regression_outcome_variables,
          function(regression_outcome_variables) {
            fit_lm(
              model = "full",
              outcome = paste(regression_outcome_variables),
              data = analysis(..1),
              tidy = TRUE
            )
          }
        )
      ),
      # add names
      lm_model_full_tidy = map(
        lm_model_full_tidy,
        ~ set_names(.,
          nm = c("PCA", "UMAP", "MEAN")
        )
      ),
      # turn list of tibbles into one tibble and create unique names in the term column required for int_pctl
      lm_model_full_tidy = map(
        lm_model_full_tidy,
        ~ map_dfr(.,
          bind_rows,
          .id = "regression_outcome"
        ) %>%
          mutate(term = paste(regression_outcome, term, sep = "."))
      )
    )
  ))


alpha_value <- c(
  "cis_90" = 0.10,
  "cis_95" = 0.05,
  "cis_99" = 0.01
)

# estimate confidence intervals
cis <- bootstraps %>%
  mutate(cis_simple = map(
    .x = boots,
    .f = ~ map(
      alpha_value,
      function(alpha_value) {
        int_pctl(
          .data = .x,
          statistics = "lm_model_simple_tidy",
          alpha = alpha_value
        )
      }
    )
  )) %>%
  mutate(cis_full = map(
    .x = boots,
    .f = ~ map(
      alpha_value,
      function(alpha_value) {
        int_pctl(
          .data = .x,
          statistics = "lm_model_full_tidy",
          alpha = alpha_value
        )
      }
    )
  )) %>%
  select(-boots, -data) %>%
  # extract names of the alpha vector
  mutate(alpha_values = map(
    .x = cis_simple,
    .f = ~ names(.x)
  )) %>%
  # unnest CIs
  unnest(c(alpha_values, cis_simple, cis_full))




# save the confidence intervals
bind_rows(
  # from the simple model
  cis %>%
    unnest(cis_simple) %>%
    select(-cis_full) %>%
    mutate(type = str_split(term, pattern = "\\.", simplify = T)[, 1]) %>%
    mutate(term = str_split(term, pattern = "\\.", simplify = T)[, 2]) %>%
    mutate(model_specification = "simple"),
  # from the full model
  cis %>%
    unnest(cis_full) %>%
    select(-cis_simple) %>%
    mutate(type = str_split(term, pattern = "\\.", simplify = T)[, 1]) %>%
    mutate(term = str_split(term, pattern = "\\.", simplify = T)[, 2]) %>%
    mutate(model_specification = "full")
) %>%
  write_csv("data/covid_fear_bootstrapped_cis.csv")

toc()
```


```{r estimation-by-country-bootstrap-confidence-intervals, eval=eval_computationally_intense}

set.seed(853147)
# for computational efficiency, we save the bootstrapped CIs in a csv
boot_times <- 5000

tic()
bootstraps_by_country <- covid_fear_nested_by_country %>%
  # sample 5000 times over each subset of the data
  mutate(boots = map(
    .x = data,
    .f = ~ rsample::bootstraps(.x,
      times = boot_times,
      apparent = T
    )
  )) %>%
  # estimate model on each of the splits
  mutate(boots = map(
    .x = boots,
    ~ mutate(
      .x,
      # estimate lm model simple
      lm_model_simple_tidy = pmap(
        list(..1 = .x$splits),
        .f = ~ map(
          regression_outcome_variables,
          function(regression_outcome_variables) {
            fit_lm(
              model = "simple",
              outcome = paste(regression_outcome_variables),
              data = analysis(..1),
              tidy = TRUE,
              by_country = TRUE
            )
          }
        )
      ),
      # add names
      lm_model_simple_tidy = map(
        lm_model_simple_tidy,
        ~ set_names(.,
          nm = c("PCA", "UMAP", "MEAN")
        )
      ),
      # turn list of tibbles into one tibble and create unique names in the term column required for int_pctl
      lm_model_simple_tidy = map(
        lm_model_simple_tidy,
        ~ map_dfr(.,
          bind_rows,
          .id = "regression_outcome"
        ) %>%
          mutate(term = paste(regression_outcome, term, sep = "."))
      ),
      # estimate lm model full
      lm_model_full_tidy = pmap(
        list(..1 = .x$splits),
        .f = ~ map(
          regression_outcome_variables,
          function(regression_outcome_variables) {
            fit_lm(
              model = "full",
              outcome = paste(regression_outcome_variables),
              data = analysis(..1),
              tidy = TRUE,
              by_country = TRUE
            )
          }
        )
      ),
      # add names
      lm_model_full_tidy = map(
        lm_model_full_tidy,
        ~ set_names(.,
          nm = c("PCA", "UMAP", "MEAN")
        )
      ),
      # turn list of tibbles into one tibble and create unique names in the term column required for int_pctl
      lm_model_full_tidy = map(
        lm_model_full_tidy,
        ~ map_dfr(.,
          bind_rows,
          .id = "regression_outcome"
        ) %>%
          mutate(term = paste(regression_outcome, term, sep = "."))
      )
    )
  ))


alpha_value <- c(
  "cis_90" = 0.10,
  "cis_95" = 0.05,
  "cis_99" = 0.01
)

# estimate confidence intervals
cis_by_country <- bootstraps_by_country %>%
  mutate(cis_simple = map(
    .x = boots,
    .f = ~ map(
      alpha_value,
      function(alpha_value) {
        int_pctl(
          .data = .x,
          statistics = "lm_model_simple_tidy",
          alpha = alpha_value
        )
      }
    )
  )) %>%
  mutate(cis_full = map(
    .x = boots,
    .f = ~ map(
      alpha_value,
      function(alpha_value) {
        int_pctl(
          .data = .x,
          statistics = "lm_model_full_tidy",
          alpha = alpha_value
        )
      }
    )
  )) %>%
  select(-boots, -data) %>%
  # extract names of the alpha vector
  mutate(alpha_values = map(
    .x = cis_simple,
    .f = ~ names(.x)
  )) %>%
  # unnest CIs
  unnest(c(alpha_values, cis_simple, cis_full))




# save the confidence intervals
bind_rows(
  # from the simple model
  cis_by_country %>%
    unnest(cis_simple) %>%
    select(-cis_full) %>%
    mutate(type = str_split(term, pattern = "\\.", simplify = T)[, 1]) %>%
    mutate(term = str_split(term, pattern = "\\.", simplify = T)[, 2]) %>%
    mutate(model_specification = "simple"),
  # from the full model
  cis_by_country %>%
    unnest(cis_full) %>%
    select(-cis_simple) %>%
    mutate(type = str_split(term, pattern = "\\.", simplify = T)[, 1]) %>%
    mutate(term = str_split(term, pattern = "\\.", simplify = T)[, 2]) %>%
    mutate(model_specification = "full")
) %>%
  write_csv("data/covid_fear_by_country_bootstrapped_cis.csv")

toc()
```


#### Model estimation

We estimate all models on the nested data. 

```{r model-estimation}

covid_fear_nested <- covid_fear_nested %>%
  mutate(
    lm_model_simple = map(
      data,
      ~ map(
        regression_outcome_variables,
        function(regression_outcome_variables) {
          fit_lm(
            model = "simple",
            outcome = paste(regression_outcome_variables),
            data = .x
          )
        }
      )
    ),
    lm_model_full = map(
      data,
      ~ map(
        regression_outcome_variables,
        function(regression_outcome_variables) {
          fit_lm(
            model = "full",
            outcome = paste(regression_outcome_variables),
            data = .x
          )
        }
      )
    )
  )



covid_fear_nested_by_country <- covid_fear_nested_by_country %>%
  mutate(
    lm_model_simple = map(
      data,
      ~ map(
        regression_outcome_variables,
        function(regression_outcome_variables) {
          fit_lm(
            model = "simple",
            by_country = TRUE,
            outcome = paste(regression_outcome_variables),
            data = .x
          )
        }
      )
    ),
    lm_model_full = map(
      data,
      ~ map(
        regression_outcome_variables,
        function(regression_outcome_variables) {
          fit_lm(
            model = "full",
            by_country = TRUE,
            outcome = paste(regression_outcome_variables),
            data = .x
          )
        }
      )
    )
  )
```

We turn the models into a tidy long dataframe with 36 observations (`covid_models`) that records each of the estimated models (6 outcome dimensions, 2 model specifications (simple vs. full), 3 dimensionality reduction methods). The `covid_models_by_country` records 72 observations because we estimate all of the models for each of the two countries separately. 

```{r regression-models-long}

covid_models <- covid_fear_nested %>%
  select(-data) %>%
  pivot_longer(
    cols = matches("model"),
    names_to = "model_specification",
    names_prefix = "lm_model_",
    values_to = "lm_model"
  ) %>%
  # model names
  mutate(lm_model = map(
    lm_model,
    ~ set_names(.,
      nm = c("PCA", "UMAP", "MEAN")
    )
  )) %>%
  # extract names of the model
  mutate(type = map(
    lm_model,
    names
  )) %>%
  unnest(c(lm_model, type))



covid_models_by_country <- covid_fear_nested_by_country %>%
  select(-data) %>%
  pivot_longer(
    cols = matches("model"),
    names_to = "model_specification",
    names_prefix = "lm_model_",
    values_to = "lm_model"
  ) %>%
  # model names
  mutate(lm_model = map(
    lm_model,
    ~ set_names(.,
      nm = c("PCA", "UMAP", "MEAN")
    )
  )) %>%
  # extract names of the model
  mutate(type = map(
    lm_model,
    names
  )) %>%
  unnest(c(lm_model, type))
```



### Figure 2: The effect of fear of COVID-19 on attitudes

We add the bootstrapped confidence intervals (CIs) to the dataframe containing the model estimates and plot the point estimates along with the respective CIs. 

```{r modelplot-estimates-covid, fig.cap="The effect of fear of COVID-19 on authoritarian, nationalist, and outgroup-hostile attitudes (left panel) and related COVID-19 policy measures (right panel). Point estimates along with 90\\%, 95\\%, and 99\\% bootstrapped percentile confidence intervals obtained from 5000 bootstrap resamples.", fig.height=6, fig.width=10}

# read bootstrapped cis
bootstrapped_cis <- read_csv("./data/covid_fear_bootstrapped_cis.csv")

scale_alpha_values <- c(
  "full" = 1,
  "simple" = 0.4
)

scale_shape_values <- c(
  "MEAN" = 17, 
  "PCA" = 15, 
  "UMAP" = 19
)

coef_plot <- covid_models %>%
  mutate(coefs = map(lm_model, tidy)) %>%
  select(-lm_model) %>%
  unnest(coefs) %>%
  select(-c("std.error", "p.value", "statistic")) %>%
  left_join(
    .,
    bootstrapped_cis %>%
      select(-.estimate)
  ) %>%
  rename(dimensionality_reduction = type) %>%
  mutate(dimensionality_reduction = str_to_upper(dimensionality_reduction)) %>%
  mutate(
    outcome =
      case_when(
        str_detect(variable, "rwa") ~ "Authoritarian",
        str_detect(variable, "nationalism") ~ "Nationalism",
        str_detect(variable, "outgroup") ~ "Outgroup-Hostile"
      ),
    level =
      case_when(
        str_detect(variable, "covid") ~ "COVID-19 policy preferences",
        TRUE ~ "Higher-level attitudes"
      )
  )



ggplot(
  coef_plot %>%
    filter(term %in% c("fear_treatment")),
  aes(
    y = outcome,
    x = estimate,
    xmin = .lower,
    xmax = .upper,
    color = fct_rev(dimensionality_reduction),
    shape = fct_rev(dimensionality_reduction), 
    alpha = model_specification,
    size = .alpha
  )
) +
  ggdist::geom_pointinterval(
    position = position_dodge(width = 0.6),
    interval_size_range = c(1, 2),
    fatten_point = 1.5
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "grey20"
  ) +
  facet_wrap(~ fct_rev(level)) +
  ylab(label = "") +
  xlab(label = "") +
  labs(caption = "") +
  jcolors::scale_color_jcolors(
    palette = "pal8",
    guide = "none", 
    name = "Method of dimensionality reduction"
  ) +
  scale_shape_manual(
    name = "Method of dimensionality reduction", 
    values = scale_shape_values
  ) +
  scale_alpha_manual(
    name = "",
    values = scale_alpha_values,
    guide = "none"
  ) + 
  labs(caption = "Note: estimates shown in light (dark) shading from simple (fully specified) models.")

```



### Regression Tables (Tables S5-S10)

In the following, we present the results shown in Figure \@ref(fig:modelplot-estimates-covid) in regression tables.


```{r modelsummary-configurations-lm-95cis}

bootstrapped95_cis <- bootstrapped_cis %>%
  filter(.alpha == .05) %>%
  nest(data = c(-variable, -model_specification, -type)) %>%
  mutate(cis = map(
    data,
    ~ tibble(
      variable = paste0(.x %>% select(term) %>% pull()),
      cis = paste0(
        "[",
        sprintf("%.2f", .x %>% pull(.lower)),
        "; ",
        sprintf("%.2f", .x %>% pull(.upper)),
        "]"
      )
    ) %>% deframe()
  )) %>%
  select(-data)


# add to the data
covid_models <- covid_models %>%
  left_join(
    .,
    bootstrapped95_cis
  ) %>%
  arrange(desc(variable), desc(model_specification))


# create regression table
cm <- c(
  "fear_treatment" = "COVID-19 Fear",
  "female1" = "Female",
  "age" = "Age",
  "urban_rural" = "Urbanity",
  "education_level" = "Education",
  "support_government" = "Gov. Support",
  "political_news" = "Pol. News",
  "religious_service" = "Church Attendance",
  "religion_life1" = "Religion Important",
  "covid_infection1" = "Covid Infection",
  "covid_vaccination1	" = "Covid Vaccination",
  "response_week_incidence" = "Covid Incidence Rate",
  "minority1" = "Minority",
  "diaspora1" = "Diaspora",
  "(Intercept)" = "Intercept"
)


# report only basic goodness-of-fit statistics
round_function <- function(x) format(round(x, 3))
gm <- list(
  list("raw" = "r.squared", "clean" = "R$^2$", "fmt" = round_function),
  list("raw" = "nobs", "clean" = "Num.Obs.", "fmt" = round_function)
)

table_caption_covid_models <- "Preferences for authoritarian, outgroup-hostile, and nationalist COVID-19 measures in response to fear of COVID-19"

table_caption_attitudes_models <- "Authoritarian, outgroup-hostile, and nationalist attitudes in response to fear of COVID-19"
```


#### PCA reduced outcome dimensions

Table \@ref(tab:pca-regression-table-fear-covid-policy) shows the effect of fear of COVID-19 on authoritarian, outgroup-hostile, and nationalist policy measures in response to the pandemic. The dependent variables are first principal components of the respective dimensions. We report confidence intervals from 5000 bootstrap resamples stratified by countries. 

Table \@ref(tab:pca-regression-table-fear-attitudes) shows the effect of fear of COVID-19 on broader levels of authoritarian, outgroup-hostile, and nationalist attitudes that are not specifically related to the pandemic. 


```{r pca-regression-table-fear-covid-policy}

table_note_covid_models_pca <- "Note: dependent variables are first principal components of the respective dimensions. 95% percentile confidence intervals from 5000 bootstrap resamples stratified by countries."


covid_models_pca <- covid_models %>%
  filter(type == "PCA") %>%
  filter(str_detect(variable,
    pattern = "covid"
  ))

names(covid_models_pca$lm_model) <-
  str_to_title(covid_models_pca$model_specification)
names(covid_models_pca$cis) <-
  str_to_title(covid_models_pca$model_specification)


modelsummary(
  models = covid_models_pca$lm_model,
  vcov = covid_models_pca$cis,
  fmt = 2,
  title = paste(
    table_caption_covid_models,
    "(outcomes: PCA)"
  ),
  coef_map = cm,
  gof_map = gm,
  escape = FALSE,
  output = "kableExtra"
) %>%
  # column labels upper dimension: pca_variable
  add_header_above(c(
    " " = 1,
    "Authoritarian" = 2,
    "Outgroup-Hostile" = 2,
    "Nationalist" = 2
  )) %>%
  # footnote
  add_footnote(table_note_covid_models_pca,
    notation = "none",
    threeparttable = TRUE
  ) %>%
  kableExtra::kable_styling(
    font_size = 10,
    latex_options = "HOLD_position"
  ) %>%
  kableExtra::row_spec(c(1, 2), background = "#eaeef0")
```


```{r pca-regression-table-fear-attitudes}

table_note_attitudes_models <- "Note: dependent variables are first principal components of the respective dimensions. 95% percentile confidence intervals from 5000 bootstrap resamples stratified by countries."


attitudes_models_pca <- covid_models %>%
  filter(type == "PCA") %>%
  filter(!str_detect(variable,
    pattern = "covid"
  ))

names(attitudes_models_pca$lm_model) <-
  str_to_title(attitudes_models_pca$model_specification)
names(attitudes_models_pca$cis) <-
  str_to_title(attitudes_models_pca$model_specification)


modelsummary(
  models = attitudes_models_pca$lm_model,
  vcov = attitudes_models_pca$cis,
  fmt = 2,
  title = paste(
    table_caption_attitudes_models,
    "(outcomes: PCA)"
  ),
  coef_map = cm,
  gof_map = gm,
  escape = FALSE,
  output = "kableExtra"
) %>%
  # column labels upper dimension: pca_variable
  add_header_above(c(
    " " = 1,
    "Authoritarian" = 2,
    "Outgroup-Hostile" = 2,
    "Nationalist" = 2
  )) %>%
  # footnote
  add_footnote(table_note_attitudes_models,
    notation = "none",
    threeparttable = TRUE
  ) %>%
  kableExtra::kable_styling(
    font_size = 10,
    latex_options = "HOLD_position"
  ) %>%
  kableExtra::row_spec(c(1, 2), background = "#eaeef0")
```


#### UMAP reduced outcome dimensions

Table \@ref(tab:umap-regression-table-fear-covid-policy) shows the effect of fear of COVID-19 on authoritarian, outgroup-hostile, and nationalist policy measures in response to the pandemic. The dependent variables are first UMAP components of the respective dimensions. We report confidence intervals from 5000 bootstrap resamples stratified by countries. 

Table \@ref(tab:umap-regression-table-fear-attitudes) shows the effect of fear of COVID-19 on broader levels of authoritarian, outgroup-hostile, and nationalist attitudes that are not specifically related to the pandemic. 



```{r umap-regression-table-fear-covid-policy}

table_note_covid_models_umap <- "Note: dependent variables are dimensions obtained by UMAP. 95% percentile confidence intervals from 5000 bootstrap resamples stratified by countries."

covid_models_umap <- covid_models %>%
  filter(type == "UMAP") %>%
  filter(str_detect(variable,
    pattern = "covid"
  ))

names(covid_models_umap$lm_model) <-
  str_to_title(covid_models_umap$model_specification)
names(covid_models_umap$cis) <-
  str_to_title(covid_models_umap$model_specification)


modelsummary(
  models = covid_models_umap$lm_model,
  vcov = covid_models_umap$cis,
  fmt = 2,
  title = paste(
    table_caption_covid_models,
    "(outcomes: UMAP)"
  ),
  coef_map = cm,
  gof_map = gm,
  escape = FALSE,
  output = "kableExtra"
) %>%
  # column labels upper dimension: umap_variable
  add_header_above(c(
    " " = 1,
    "Authoritarian" = 2,
    "Outgroup-Hostile" = 2,
    "Nationalist" = 2
  )) %>%
  # footnote
  add_footnote(table_note_covid_models_umap,
    notation = "none",
    threeparttable = TRUE
  ) %>%
  kableExtra::kable_styling(
    font_size = 10,
    latex_options = "HOLD_position"
  ) %>%
  kableExtra::row_spec(c(1, 2), background = "#eaeef0")
```


```{r umap-regression-table-fear-attitudes}

table_note_attitudes_models_umap <- "Note: dependent variables are dimensions obtained by UMAP. 95% percentile confidence intervals from 5000 bootstrap resamples stratified by countries."


attitudes_models_umap <- covid_models %>%
  filter(type == "UMAP") %>%
  filter(!str_detect(variable,
    pattern = "covid"
  ))

names(attitudes_models_umap$lm_model) <-
  str_to_title(attitudes_models_umap$model_specification)
names(attitudes_models_umap$cis) <-
  str_to_title(attitudes_models_umap$model_specification)


modelsummary(
  models = attitudes_models_umap$lm_model,
  vcov = attitudes_models_umap$cis,
  fmt = 2,
  title = paste(
    table_caption_attitudes_models,
    "(outcomes: UMAP)"
  ),
  coef_map = cm,
  gof_map = gm,
  escape = FALSE,
  output = "kableExtra"
) %>%
  # column labels upper dimension: umap_variable
  add_header_above(c(
    " " = 1,
    "Authoritarian" = 2,
    "Outgroup-Hostile" = 2,
    "Nationalist" = 2
  )) %>%
  # footnote
  add_footnote(table_note_attitudes_models_umap,
    notation = "none",
    threeparttable = TRUE
  ) %>%
  kableExtra::kable_styling(
    font_size = 10,
    latex_options = "HOLD_position"
  ) %>%
  kableExtra::row_spec(c(1, 2), background = "#eaeef0")
```


#### Means reduced outcome dimensions

Table \@ref(tab:mean-regression-table-fear-covid-policy) shows the effect of fear of COVID-19 on authoritarian, outgroup-hostile, and nationalist policy measures in response to the pandemic. The dependent variables are simple means of all items belonging to the respective dimensions. We report confidence intervals from 5000 bootstrap resamples stratified by countries. 

Table \@ref(tab:mean-regression-table-fear-attitudes) shows the effect of fear of COVID-19 on broader levels of authoritarian, outgroup-hostile, and nationalist attitudes that are not specifically related to the pandemic. 


```{r mean-regression-table-fear-covid-policy}

table_note_covid_models_mean <- "Note: dependent variables are dimensions obtained by mean aggregation. 95% percentile confidence intervals from 5000 bootstrap resamples stratified by countries."

covid_models_mean <- covid_models %>%
  filter(type == "MEAN") %>%
  filter(str_detect(variable,
    pattern = "covid"
  ))

names(covid_models_mean$lm_model) <-
  str_to_title(covid_models_mean$model_specification)
names(covid_models_mean$cis) <-
  str_to_title(covid_models_mean$model_specification)


modelsummary(
  models = covid_models_mean$lm_model,
  vcov = covid_models_mean$cis,
  fmt = 2,
  title = paste(
    table_caption_covid_models,
    "(outcomes: mean)"
  ),
  coef_map = cm,
  gof_map = gm,
  escape = FALSE,
  output = "kableExtra"
) %>%
  # column labels upper dimension: mean_variable
  add_header_above(c(
    " " = 1,
    "Authoritarian" = 2,
    "Outgroup-Hostile" = 2,
    "Nationalist" = 2
  )) %>%
  # footnote
  add_footnote(table_note_covid_models_mean,
    notation = "none",
    threeparttable = TRUE
  ) %>%
  kableExtra::kable_styling(
    font_size = 10,
    latex_options = "HOLD_position"
  ) %>%
  kableExtra::row_spec(c(1, 2), background = "#eaeef0")
```


```{r mean-regression-table-fear-attitudes}

table_note_attitudes_models_mean <- "Note: dependent variables are dimensions obtained by mean aggregation. 95% percentile confidence intervals from 5000 bootstrap resamples stratified by countries."


attitudes_models_mean <- covid_models %>%
  filter(type == "MEAN") %>%
  filter(!str_detect(variable,
    pattern = "covid"
  ))

names(attitudes_models_mean$lm_model) <-
  str_to_title(attitudes_models_mean$model_specification)
names(attitudes_models_mean$cis) <-
  str_to_title(attitudes_models_mean$model_specification)


modelsummary(
  models = attitudes_models_mean$lm_model,
  vcov = attitudes_models_mean$cis,
  fmt = 2,
  title = paste(
    table_caption_attitudes_models,
    "(outcomes: mean)"
  ),
  coef_map = cm,
  gof_map = gm,
  escape = FALSE,
  output = "kableExtra"
) %>%
  # column labels upper dimension: mean_variable
  add_header_above(c(
    " " = 1,
    "Authoritarian" = 2,
    "Outgroup-Hostile" = 2,
    "Nationalist" = 2
  )) %>%
  # footnote
  add_footnote(table_note_attitudes_models_mean,
    notation = "none",
    threeparttable = TRUE
  ) %>%
  kableExtra::kable_styling(
    font_size = 10,
    latex_options = "HOLD_position"
  ) %>%
  kableExtra::row_spec(c(1, 2), background = "#eaeef0")
```

### Country-specific results

Figure \@ref(fig:modelplot-estimates-covid-by-country) presents the estimated effects of COVID-19 related anxieties on the different outcome dimensions which have been obtained by fitting the 36 regressions shown in \@ref(fig:modelplot-estimates-covid) separately for each of the two countries.

```{r modelplot-estimates-covid-by-country, fig.cap="The effect of fear of COVID-19 on authoritarian, nationalist, and outgroup-hostile attitudes (left panel) and related COVID-19 policy measures (right panel). Results obtained by country-specific regressions. Point estimates along with 90\\%, 95\\%, and 99\\% bootstrapped percentile confidence intervals obtained from 5000 bootstrap resamples.", fig.height=8, fig.width = 10}

bootstrapped_cis_by_country <- read_csv("./data/covid_fear_by_country_bootstrapped_cis.csv")

scale_alpha_values <- c(
  "full" = 1,
  "simple" = 0.4
)

scale_shape_values <- c(
  "MEAN" = 17, 
  "PCA" = 15, 
  "UMAP" = 19
)

coef_plot_by_country <- covid_models_by_country %>%
  mutate(coefs = map(lm_model, tidy)) %>%
  select(-lm_model) %>%
  unnest(coefs) %>%
  select(-c("std.error", "p.value", "statistic")) %>%
  left_join(
    .,
    bootstrapped_cis_by_country %>%
      select(-.estimate)
  ) %>%
  rename(dimensionality_reduction = type) %>%
  mutate(dimensionality_reduction = str_to_upper(dimensionality_reduction)) %>%
  mutate(
    outcome =
      case_when(
        str_detect(variable, "rwa") ~ "Authoritarian",
        str_detect(variable, "nationalism") ~ "Nationalism",
        str_detect(variable, "outgroup") ~ "Outgroup-Hostile"
      ),
    level =
      case_when(
        str_detect(variable, "covid") ~ "COVID-19 policy preferences",
        TRUE ~ "Higher-level attitudes"
      )
  )



ggplot(
  coef_plot_by_country %>%
    filter(term %in% c("fear_treatment")),
  aes(
    y = outcome,
    x = estimate,
    xmin = .lower,
    xmax = .upper,
    color = fct_rev(dimensionality_reduction),
    shape = fct_rev(dimensionality_reduction), 
    alpha = model_specification,
    size = .alpha
  )
) +
  ggdist::geom_pointinterval(
    position = position_dodge(width = 0.6),
    interval_size_range = c(1, 2),
    fatten_point = 1
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "grey20"
  ) +
  facet_grid(
    cols = vars(fct_rev(level)),
    rows = vars(country)
  ) +
  ylab(label = "") +
  xlab(label = "") +
  labs(caption = "") +
  jcolors::scale_color_jcolors(
    palette = "pal8",
    name = "Method of dimensionality reduction", 
    guide = "none"
  ) +
  scale_alpha_manual(
    name = "",
    values = scale_alpha_values,
    guide = "none"
  ) + 
  scale_shape_manual(
    name = "Method of dimensionality reduction", 
    values = scale_shape_values
  )  +
  labs(caption = "Note: estimates shown in light (dark) shading from simple (fully specified) models.")

```


### Within-dimensions analyses

We also report the simple differences in means between treatment and control across the various items of the different dimensions. Figure \@ref(fig:within-country-within-dimensions) shows that there are no statistically significant differences on any of the outcome items (for a detailed description of the items see Tables SI1 and SI2 in the Supplementary Material), neither for the Hungarian respondents, nor for the Romanian respondents. 

The graph shows that the mean values on the respective outcome variables (standardised to a zero mean and unit standard variation) among those respondents to whom their fears and anxieties during the COVID-19 pandemic were cognitively accessible ("Treatment") are statistically indistinguishable from the mean values among those respondents to whom their fears and anxieties during the COVID-19 pandemic were not cognitively accessible ("Control"). 

```{r within-country-within-dimensions, fig.height = 7, fig.width=9, fig.cap="The effect of fear recall on the various outcome items within each dimension."}

# create a long dataset with all the different outcomes
item_outcomes_long <- covid_fear_data %>%
  select(
    country, fear_treatment, ID,
    all_of(outcome_variables)
  ) %>%
  group_by(country) %>%
  mutate(across(all_of(outcome_variables), ~ as.numeric(scale(.x)))) %>%
  pivot_longer(
    cols = -c("ID", "country", "fear_treatment"),
    names_to = "outcome_variable",
    values_to = "value"
  ) %>%
  # rename the variables for the graph
  mutate(outcome_group = case_when(
    str_detect(outcome_variable, "covid_policy_outgroup_") ~ "outgroup_covid",
    str_detect(outcome_variable, "covid_policy_nationalism_") ~ "nationalism_covid",
    str_detect(outcome_variable, "covid_policy_authoritarianism") ~ "rwa_covid",
    str_detect(outcome_variable, "rwa_") ~ "rwa",
    str_detect(outcome_variable, "outgroup_") ~ "outgroup",
    str_detect(outcome_variable, "nationalism_") ~ "nationalism"
  )) %>%
  mutate(level = case_when(
    str_detect(outcome_variable, "covid") ~ "COVID-19 policy preferences",
    TRUE ~ "Higher-level attitudes"
  )) %>%
  # rename the variables
  mutate(outcome_variable = case_when(
    str_detect(outcome_variable, "rwa") ~ str_replace(
      outcome_variable,
      "rwa", "Authoritarian"
    ),
    str_detect(outcome_variable, "covid_policy_") ~ str_remove_all(
      outcome_variable,
      "covid_policy_"
    ),
    TRUE ~ outcome_variable
  )) %>%
  mutate(outcome_variable = case_when(
    str_detect(outcome_variable, "nationalism") ~ str_replace(
      outcome_variable,
      "nationalism",
      "Nationalist"
    ),
    str_detect(outcome_variable, "outgroup") ~ str_replace(
      outcome_variable,
      "outgroup",
      "Outgroup-Hostile"
    ),
    str_detect(outcome_variable, "authoritarianism") ~ str_replace(
      outcome_variable,
      "authoritarianism",
      "Authoritarian"
    ),
    TRUE ~ outcome_variable
  )) %>%
  mutate(fear_treatment = ifelse(fear_treatment == 1,
    "Treatment",
    "Control"
  ))


ggplot(
  data = item_outcomes_long,
  aes(
    x = value,
    y = outcome_variable,
    color = as.factor(fear_treatment)
  )
) +
  facet_grid(
    cols = vars(country),
    rows = vars(level),
    scales = "free_y"
  ) +  
  stat_summary(
    fun.data = "mean_cl_boot",
    B = 5000,
    geom = "pointrange",
    position = position_dodge(width = 0.5)
  ) +
  xlab("") +
  ylab("") +
  jcolors::scale_color_jcolors(
    name = "",
    palette = "pal9"
  )
```

# Session Info

- `r paste0(benchmarkme::get_sys_details()$sys_info$sysname, " (", benchmarkme::get_sys_details()$sys_info$version, ")")` 
- `r paste0(benchmarkme::get_cpu()$model_name, " with ", benchmarkme::get_cpu()$no_of_cores,  " cores")`
- `r R.Version()$version.string`


```{r session-info}

xfun::session_info()

```


